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APPROXIMATELY UNBIASED TESTS OF REGIONS USING 
MULTISTEP-MULTISCALE BOOTSTRAP RESAMPLING^ 

By Hidetoshi Shimodaira 
Tokyo Institute of Technology 

Approximately unbiased tests based on bootstrap probabilities 
are considered for the exponential family of distributions with un- 
known expectation parameter vector, where the null hypothesis is 
represented as an arbitrary-shaped region with smooth boundaries. 
This problem has been discussed previously in Efron and Tibshirani 
[Ann. Statist. 26 (1998) 1687-1718], and a corrected p- value with 
second-order asymptotic accuracy is calculated by the two-level boot- 
strap of Efron, Halloran and Holmes [Proc. Natl. Acad. Sci. U.S.A. 
93 (1996) 13429-13434] based on the ABC bias correction of Efron 
[J. Amer. Statist. Assoc. 82 (1987) 171-185]. Our argument is an 
extension of their asymptotic theory, where the geometry, such as 
the signed distance and the curvature of the boundary, plays an im- 
portant role. We give another calculation of the corrected p-value 
without finding the "nearest point" on the boundary to the obser- 
vation, which is required in the two-level bootstrap and is an im- 
plementational burden in complicated problems. The key idea is to 
alter the sample size of the replicated dataset from that of the ob- 
served dataset. The frequency of the replicates falling in the region 
is counted for several sample sizes, and then the p-value is calcu- 
lated by looking at the change in the frequencies along the changing 
sample sizes. This is the multiscale bootstrap of Shimodaira [Sys- 
tematic Biology 51 (2002) 492-508], which is third-order accurate for 
the multivariate normal model. Here we introduce a newly devised 
multistep- multiscale bootstrap, calculating a third-order accurate p- 
value for the exponential family of distributions. In fact, our p- value is 
asymptotically equivalent to those obtained by the double bootstrap 
of Hall [The Bootstrap and Edge-worth Expansion (1992) Springer, 
New York] and the modified signed likelihood ratio of Barndorff- 
Nielsen [Biometnka 73 (1986) 307-322] ignoring 0(n"^/^) terms, yet 
the computation is less demanding and free from model specification. 
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The algorithm is remarkably simple despite complexity of the theory 
behind it. The differences of the p- values are illustrated in simple ex- 
amples, and the accuracies of the bootstrap methods are shown in a 
systematic way. 

1. Introduction. We start with a simple example of Efron and Tibshirani 
(1998) to illustrate the issue to discuss. Let Xi, . . . , Xn be independent p- 
dimensional multivariate normal vectors with mean vector and covariance 
matrix identity Ip, 

Xi, . . . ,X„ ~ Np{fi, Ip). 

For given observed values let us assume that we would like to 

know whether = /i^ + • • • + /i^ < 1 or not. The problem is also de- 
scribed in a transformed variable Y = ^/nX with mean r/ = \/n[i^ where x = 
{x\ + • • • + Xn)/n is the sample average. We have observed a p-dimensional 
multivariate normal vector y having unknown mean vector rj and covariance 
matrix the identity, 

(1.1) Y^Npi^Jp). 

Then the null hypothesis we are going to test is rj gTZ, with the spherical 
region 

(1.2) n = {7]:\\7]\\<^}. 

This problem is simple enough to give the exact answer. The frequentist 
confidence level, namely, the probability value (p-value) for the spherical 
null hypothesis is calculated as the probability of \\Y\\'^ being greater than 
or equal to the observed assuming that t] is on the boundary dTZ = 
{77: \\r]\\ = \/n} of TZ. The exact p-value is easily calculated knowing that 
is distributed as the chi-square distribution with degrees of freedom p 
and noncentrality 

In this paper we are going to remove two restrictions in the above problem 
for generalization, (i) The underlying probability model for Y is the expo- 
nential family of distributions, instead of the multivariate normal model; we 
denote the density function with the expectation parameter rj as 

(1.3) Y^f{y;ri). 

(ii) The null hypothesis will be represented as an arbitrarily-shaped region 
TZ with smooth boundaries, instead of the spherical region. The surface of 
dTZ may be represented as the Taylor series with coefhcients d"^, e"'^'^, . . . 

(1.4) A7]p = -d'^'AT^aArib - e'''''ArjaArjbAri, + ■■■ 

in the local coordinates (Arji, . . . , Arjp) by taking the origin at a point on 
dTZ and rotating the axes properly. The summation convention such as 
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d^'' ArjaA-Vb = J2a=i Efe=i d'^^ ^Va^Tlb will be used, where the indices a,b,. . . 
may run through 1, ... ,p — 1 and . . . may run though 1,. . . ,p when used 
as subscripts or superscripts for p-dimensional vectors. The axes are taken 
so that Ar/i, . . . , Arjp-i are for the tangent space of the surface, and Aijp is 
for its orthogonal space taken positive in the direction pointing away from 
TZ. This general setting is the "problem of regions" discussed previously 
in Efron and Tibshirani (1998), and our argument is an extension of their 
asymptotic theory, where the geometry, such as the signed distance and the 
curvature of the boundary, plays an important role. 

Since the exact p-value is available only for special cases, we will discuss 
several bootstrap methods to calculate approximate p-values from y under 
the assumptions (i) and (ii) above. Let a denote a specified significance 
level, and a{y) denote an approximate p-value. A large value of a{y) may 
indicate evidence to support the null hypothesis rj gTZ. On the other hand, 
if a(y) < a is observed, then we reject the null hypothesis and conclude 
that rj ^TZ. The hypothesis test of TZ is said to be unbiased if the rejection 
probability is equal to a whenever rj £ dTZ. The approximate p- value is said 
to be kth order accurate if the asymptotic bias is of order 0{n~^/'^), that is, 

(1.5) Pr{d(y) <a;r/} = a + 0(n~'=/2)^ r/ G 57^, 

holds for < a < 1. For sufficiently large n, approximately unbiased p- values 
of higher-order accuracy are considered to be better than those of lower-order 
accuracy. 

We will not specify the probabilistic model or the shape of the region ex- 
plicitly in the calculation of the p-value, but only assume that a mechanism 
is available to us for generating the bootstrap replicates and identifying 
whether the outcomes are in the region or not. This setting is important 
for complicated practical applications, where the exact p-value is not avail- 
able and, thus, bootstrap methods are used for approximation. The phylo- 
genetic tree selection discussed in Efron, Halloran and Holmes (1996) and 
Shimodaira (2002) is a typical case; the history of evolution represented as a 
tree is inferred by a model-based clustering of the DNA sequences of organ- 
isms, where we are given complex computer software for inferring the tree 
from a dataset. For calculating p- values of the hypothetical evolutionary 
trees, we can easily run bootstrap simulations, although computationally 
demanding, by repeatedly applying the software to replicated datasets. 

We confine our attention to the parametric bootstrap of continuous ran- 
dom vectors for mathematical simplicity. We also assume that the boundary 
of the region is a smooth surface. In practical applications, however, it is 
often the case that the nonparametric bootstrap is employed, the random 
vector is discrete and the boundary is nonsmooth. Regions with nonsmooth 
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boundaries, in particular, may lead to serious difficulty as discussed in Perl- 
man and Wu (1999, 2003). Further study is needed to bridge these gaps 
between the theory and practice. 

The frequency of the bootstrap replicates falling in the region, namely, 
the bootstrap probability, has been used widely since its application to phy- 
logenetic tree selection in Felsenstein (1985). This is also named "empirical 
strength probability" of TZ in Liu and Singh (1997), where a modification for 
nonsmooth boundary is discussed as well. The bootstrap probability is, how- 
ever, biased as an approximation to the exact p- value and, thus, the two-level 
bootstrap of Efron, Halloran and Holmes (1996) and Efron and Tibshirani 
(1998) is developed to improve the accuracy. Under the assumptions (i) 
and (ii) above, the two-level bootstrap calculates a second-order accurate 
p-value, whereas the bootstrap probability is only first-order accurate. 

The bias of the bootstrap probability mainly arises from the curvature 
of dTZ. The two-level bootstrap estimates the curvature for bias correc- 
tion, where the curvature is estimated by generating second-level replicates 
around fi(y). Here fi{y) denotes the maximum likelihood estimate for rj re- 
stricted to dTZ. fj{y) is the nearest point on dTZ to y for (1.1). For the spher- 
ical region, fi{y) = ^/ny/\\y\\ is easily obtained, but 7]{y) must be obtained 
by numerical search in general, leading to an implementational burden in 
complex problems. This motivated our development of a new method. 

The multiscale bootstrap is developed in Shimodaira (2002) to calculate 
another bias corrected p- value. It does not require fj{y). Instead, the boot- 
strap probabilities are calculated for sets of bootstrap replicates with several 
sample sizes which may differ from that of the observed data. This, in effect, 
alters the scale parameter of the replicates (Figure 1). The key idea is to 
estimate the curvature from the change in the bootstrap probabilities along 
varying sample sizes. The corrected p-value is third-order accurate for any 
arbitrarily-shaped region with smooth boundaries under the multivariate 
normal model. The normality assumption is not as restrictive as it might 
look at first, because the procedure is transformation-invariant and should 
work fine if there exists a transformation from the dataset to the normal Y 
and if the null hypothesis is represented as a region of r/. We do not have 
to know what the transformation is. However, it becomes only first-order 
accurate if there is no such transformation to (1.1) but only one to (1.3). 

The multiscale bootstrap can be used easily for complex problems. It is as 
easy as the usual bootstrap. We only have to change the sample size of the 
bootstrap replicates, and apply a regression fit to the bootstrap probabilities. 
The bias corrected p- value is calculated from the slope of the regression curve 
(Figure 2). This procedure is implemented in computer software [Shimodaira 
and Hasegawa (2001)] for phylogenetic tree selection, and is also applied to 
gene network estimation from microarray expression profiles [Kamimura et 
al. (2003)]. In these applications, the multiscale bootstrap can calculate the 
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p-values for many related hypotheses at the same time; we do not have to 
run time-consuming bootstrap simulations separately for these hypotheses. 
For example, biologists are interested in the monophyletic hypothesis that 
some specified species constitute a cluster in the phylogenetic tree, and there 
are many such hypotheses for groups of species. The bootstrap probabilities 
for these hypotheses are obtained at the same time from a single run of 
bootstrap simulation for each scale. We only have to apply the regression fit 
separately to the multiscale bootstrap probabilities of each hypothesis. 

In this paper we provide the theoretical foundation of the multiscale boot- 
strap, and introduce a newly devised multistep-multiscale bootstrap resam- 
pling. This method calculates an approximately unbiased p- value with third- 
order asymptotic accuracy under the assumptions (i) and (ii) . The previously 
developed method of Shimodaira (2002) corresponds to a special case of the 
new method, that is, the one-step multiscale bootstrap. 

For explaining the bootstrap methods, a rather intuitive argument is given 
in Sections 2 to 6 using simple examples. A more formal argument is given 
in Section 7, and the technical details are given in a supporting document 
[Shimodaira (2004)]. We introduce a modified signed distance, and give a 
unified approach to the asymptotic analysis of the bootstrap methods using 
Edgeworth series, as well as the tube formula of Weyl (1939). Third-order 
accuracy is also shown there for the p- value computed by the modified signed 
likelihood ratio [Barndorff'-Nielsen (1986)], which requires the analytic ex- 
pression of the likelihood function, and for the p-value computed by the 




Fig. 1. Multiscale bootstrap. The three circles with dashed lines indicate the conditional 
distributions of the bootstrap replicates with mean y and scales r = l/\/2, 1, V^- In this 
particular configuration, the bootstrap probability may increase by halving the sample size 
to alter t — 1 to \f2, and may decrease by doubling the sample size to alter t = 1 to 1/ \f2. 
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double bootstrap [Hall (1992)], which requires a huge number of replicates, 
as well as computation of ?)(y). The multistep-multiscale bootstrap method 
requires only the bootstrap mechanism for generating replicates around y, 
inheriting the simplicity from the one-step multiscale bootstrap. The price 
for higher-order accuracy and simpler implementation is a large number of 
replicates, which can be as large as that of the double bootstrap. These three 
p-values are, in fact, shown to be equivalent ignoring terms. 

Our argument may not be justified unless the assumptions (i) and (ii) 
hold. We are not sure yet how robust the multistep-multiscale bootstrap 
method is under misspecifications of the exponential family model. It is 
shown at the end of Section 4, however, that the one-step method adjusts 
the bias halfway, though not completely, under misspecifications of the nor- 
mal model. A simulation study in Shimodaira (2002) shows that the bias of 
the one-step method under the normal model is very small even if the bound- 
ary is piecewise smooth, but the bias becomes larger as rj moves closer to 
nonsmooth points on the boundary. 

2. Two-level bootstrap resampling. Although our ultimate goal is to get 
rid of the normal assumption, we use normality in this section to illustrate 
the bootstrap methods, and besides (1.1), we also assume (1.2). For given 
observed value x, we consider the parametric bootstrap resampling 



Typically, the sample size ni of the replicated dataset should be equal to n, 
but we reserve the generality of using any value for ni . The scaling factor of 
the bootstrap, ri = ^/n/ni, will be altered later in the multiscale bootstrap. 
Once we specify ti, we may generate B, say 10,000, replicated datasets, and 
compute the average X = {Xi + • • • + X*J/ni for each replicate. A large 
value of the frequency that H-'^*!!^ < 1 holds in the replicates may indicate 
a high chance of the null hypothesis < 1 being correct. This is also 
described in a transformed variable Y* = ^/nX . For given observed value 
y, we consider the parametric bootstrap resampling 



where the index 1 indicates the "one-step" bootstrap in connection with a.^ 
and defined later, as shown in Table 1. a\ is estimated by the frequency 
of y* G 7?. from the B bootstrap replicates with the binomial variance o.\{\ — 



Xl,...,X:^r^Np{x,Ip). 




Y* ^Np{y,Tflp), 



and the bootstrap probability with scale ri is denoted by 



dl(y,rl)=Pr{y*G7^;2/,rl} 



ai)/B. 
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Let us consider a numerical example with 

(2.2) p = A, n=W, ||xf = 2.680. 

Although 1 1 2; IP > 1, we are not sure if < 1 holds or not. The frequentist 
confidence level for the null hypothesis is given by the exact p- value, which 
we will denote by aoo{y), or simply a^o for brevity sake. In this numerical 
example, the value of is, in fact, chosen to make aoo{y) = 0.05. cioo may 
be approximated by the bootstrap probability with ri = 1, denoted by 

ao(y) = ai(y,l). 

This turns out to be ao(y) = 0.0085, showing qq is not a very good approx- 
imation to doo- Here the problem is so simple that q;o(?/)) as well as aoo{y), 
can be computed numerically from the noncentral chi-square distribution 
function. If the bootstrap resampling with B = 10,000, say, is used for 6lq, 
the standard error becomes 0.0009. 

A modification of do is developed based on the geometric theory in Efron, 
Halloran and Holmes (1996) and Efron and Tibshirani (1998) to improve the 
accuracy of the approximation to doo- The idea is to compute do(?7(y)) by 
generating the second-level replicates around fi{y) for estimating the cur- 
vature of the surface dTZ. When the surface of dTZ is flat, do(??(y)) = \. 
It becomes smaller/larger than ^ when the surface is curved toward/away 
from TZ. Let z denote a generic symbol for the z-value corresponding to a 
p-value a with relation z = — $~^(a), where $~^(-) is the inverse of the 
standard normal distribution function $(•). For example, we may write 
-^o(y) = — ^~^("o(y))- The ABC conversion formula of Efron (1987) and 
DiCiccio and Efron (1992) is 

(2-3) 2abc(y) = r. ... , . — ^TTTTY - zoiviy)), 

Table 1 

Bootstrap probabilities and corrected p-values 



Symbol 


Section 


Description 


&i{y,Ti) 


2 


Bootstrap probability 


&oo{y) 


2 


Exact p-value* 


«o(y) 


2 


Bootstrap probability (ri = 1) 




2 


Two-level bootstrap corrected p-value 


aiiy) 


3 


Multiscale bootstrap corrected p-value 


0!2{y,ri,T2) 


4 


Two-step bootstrap probability 




4 


Two-step multiscale bootstrap corrected p-value 


0!3iy,ri,T2,T3) 


5 


Three-step bootstrap probability 


asiy) 


5 


Three-step multiscale bootstrap corrected p-value 



*A third-order accurate p-value in Section 7. 
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where 5abc(y)i zo{y), and zo{f]{y)) are denoted Z, Z, and respectively, in 
the notation of equation (6.6) of Efron and Tibshirani (1998). The corrected 
p- value for the two-level bootstrap is then defined by aabc(?/) = ^(—Zabciy))- 
The acceleration constant a, characterizing the probabilistic model, is known 
to be a = for the normal model, a may also be estimated using the second- 
level bootstrap for (1.3); for details we refer to Efron, Halloran and Holmes 
(1996). Note that the sign in front of a in (2.3) is reversed from that of 
equation (6.6) of Efron and Tibshirani (1998), because the Ar/p-axis is taking 
the opposite direction here. 

The p-values for the numerical example of (2.2) are 

Qo(y) = 0.0085, ao(??(y)) = 0.315, 

aabc(y) = 0.0775, aoo(y) = 0.05. 

We observe that Aabc shows great improvement over uq to approximate ctoo • 
This improvement is also confirmed in the asymptotic argument. It has been 
shown in Efron and Tibshirani (1998) that = 1 for ao, and k = 2 for Oabc 
under (1.3) and (1.4). 

3. Multiscale bootstrap resampling. Here we continue to use the normal 
model (1.1) for the argument of the corrected p-value in this section. The 
bootstrap probability changes if the replicate sample size changes. When we 
alter ni = 10 to ni = 3 for the numerical example of (2.2), or equivalently al- 
ter the scale ti = 1 to ri = ^10/3, we observe that ai{y, 1) = 0.0085 changes 
to ai{y, \/W/3) = 0.0359. In the multiscale bootstrap, ai(y, ri) is computed 
for several values of ri = \Jnjn\. For example, instead of n = 10, we use the 
following five n\ values: 

(3.1) ni =3,6,10,15,21, 

and compute the corresponding bootstrap probabilities 

(3.2) ai(y, n) = 0.0359, 0.0205, 0.0085, 0.0028, 0.0008. 

These values, as well as those for other parameter settings, are shown in 
Figure 2 by plotting the z-value along the inverse of the scale. The hori- 
zontal axis is 1/ri = ^Jni/n = 0.55, 0.78, 1, 1.23, 1.45, and the vertical axis is 
5i(2/,Ti) = -$-i(ai(y, n)) = 1.80, 2.04, 2.39, 2.77, 3.17. 

Figure 2 shows these values along with a regression fit. This is obtained 
by fitting a regression model with explanatory variables 1/ri and ri, 

(3.3) zi{y,Ti)^v/Ti+CTi, 

to the plot, where v and c are the regression coefficients estimated as 

(3.4) i) = 2.002, c = 0.385 

for the plot of (3.2). We observe that the regression fit agrees with the plots 
very well for the cases in Figure 2. The regression model (3.3) has been 
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justified in Shimodaira (2002) under (1.1) and (1.4); we will use "~" to 
indicate that equality holds up to 0(n~^) terms with the error of order 
0(n-3/2). The regression model with explanatory variables 1/ri and ri will 
be justified later, in fact, under (1.3) and (1.4) as seen in (7.15), although the 
following interpretation of the coefficients should be modified accordingly. 

A simple geometric interpretation can be given to the regression coef- 
ficients under (1.1) and (1.4). Efron and Tibshirani (1998) have shown a 
formula equivalent to 

(3.5) zo{y)^v + c, 

where v and c correspond to xq and di — xod2, respectively, in their equa- 
tion (2.19). V is the signed distance of Efron (1985), defined as the distance 
from y to dTZ with a positive/negative sign when y is outside/inside of TZ. 
Thus, V = ib||y — 77(?/)|| measures evidence of the null hypothesis being wrong, 
c is related to the (p — 1) x (p — 1) matrix d""^ measuring the curvature of 
dlZ at ?}(?/); d"'^ is defined as d""^ in (1.4) by making the local coordinates 
orthonormal at r}(y). In our notation, c = di — vd2, where di = d"''^ is the 
trace of and = {d^^f = Ea=l iXXid''^? is that for the squared ma- 
trix. When dlZ is flat at 77(7/), d""^ = and, thus, c = 0. v, di and ^2 are 
transformation-invariant functions of y calculated from the shape of the 
boundary and the density function of Y; they are referred to as geometric 



a« =0.05 




0.6 a.s 1 1.2 1.4 
1/r 



a«. =0.95 




J I I I L 



0,6 0.8 1 1,2 1.4 

1/T 



Fig. 2. Plots of the z-value of the multiscale bootstrap probability along the inverse of the 
scale T for the normal example (p = 4) of Section 2 and the exponential example (p — 1) of 
Section 4. Parameter values are chosen so that the exact p-value is either 0.05 (left panel) 
or 0.95 (right panel). The curves are drawn by the regression model of equation (3.3). 
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quantities here. Under (1.1) and (1.2) these quantities are 

(3.6) v = \\y\\-^, d^ = E^l^ d^ = P^_ 

This computes directly, 

(3.7) ?) = 2.015, c = 0.323 

for (2.2), showing good agreement with those computed indirectly from the 
multiscale bootstrap, v and c in (3.4) are actually estimating those in (3.7), 
thus, it would be appropriate to denote the former as v and c, although we do 
not make the notational distinction. This estimation is third-order accurate, 
since the regression model (3.3) holds for (3.7) with error of 0(n~3/2). 
Considering that v and c are functions of y, we may define a statistic 

(3.8) h{.y)=v-c. 

This is equivalent to the pivot statistic of Efron (1985), and Pr{zi(y) < 
x\r]\ ~ ^(x) for T] G dTt under (1.1) and (1.4); see equation (2.16) of Efron 
and Tibshirani (1998). Thus, a third-order accurate value is defined by 
di(y) = $(— zi(y)). We can compute ai(y) using v and c obtained from the 
multiscale bootstrap. For the example of (2.2), 

ai(y) = $(-2,002 0.385) = 0.0529, 

showing an improvement over aabc(y) = 0.0775 to approximate aoo(y) = 
0.05. The index of a\ indicates the "one-step" bootstrap as similarly for a\. 

It is interesting to note that we can also read off the values of z\ (y) from 
Figure 2. The differentiation of (3.3) with respect to 1/ri is 

dzi{y,Ti) , ^ 2 
9(l/ri) 

and the slope of the regression curve at 1/ti = 1 gives z\{y). The corrected 
p-value a\ is essentially obtained from the change of the bootstrap proba- 
bility in the multiscale bootstrap. 



4. Two-step multiscale bootstrap resampling. The one-step multiscale 
bootstrap described in Section 3 calculates a very accurate p-value for the 
arbitrarily-shaped region if there exists a transformation from the dataset 
to the normal model. However, it can be inaccurate if such a transformation 
does not exist even approximately. This restriction essentially comes from 
the fact that the covariance matrix of y in (1.1) is constant with respect 
to T]. The acceleration constant a of the ABC formula measures the rate of 
change in the covariance matrix, and a is assumed zero in the derivation of 
(3.8). Here we introduce the two-step multiscale bootstrap for estimating a 
to improve the accuracy of the one-step multiscale bootstrap. 
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A breakdown of the one-step multiscale bootstrap method is illustrated 
in the following example. Let Xi, . . . ,Xn be one-dimensional independent 
exponential random variables with mean /i, 

Xi, . . . ,X„ ~ exp(-x//i - log/i), 

and let the null hypothesis of interest be /U < 1. The exact p- value is cal- 
culated by knowing that a transformed variable Y = ^JnX is distributed as 
Gamma with shape n and mean r] = -y/n^. We consider a numerical example 
with 

(4.1) p = l, n = 10, x = 1.571, 

so that ftoo (y) = 0.05. The multiscale bootstrap probabilities for the five n\ 
values in (3.1) are computed as 

(4.2) ai(y, n) = 0.2990, 0.1875, 0.1115, 0.0622, 0.0322, 

and the regression coefficients of (3.3) are estimated as u = 1.328, c = —0.110. 
Then the corrected p-value is computed as 

(4.3) di(y) = ^>(-1.328 - 0.110) = 0.0753. 

Although this is an improvement over a;o(?/) = 0.112, it is not as good as 
in the normal example above. The pivot (3.8) is not justified under (1.3) in 
general, and t3;i(7/) is, in fact, only first-order accurate for the exponential 
example. 

The two-step multiscale bootstrap is employed simply to generate a second- 
step replicate from every first-step replicate. Let us denote the conditional 
density of the first-step bootstrap replicate Y* = ^/nX as 

(4.4) y*~/(y*;y,ri), 

given mean y = yJnX and scale t\ under (1.3), which reduces to f{y*;y, 1) = 
f{y*;y) when ri = yjnjnx is unity. This becomes (2.1) for (1.1), and Gamma 
with shape n\ and mean y for the exponential example. We generate a 
second-step replicate Y** for each y* . The conditional density of Y** given 
y* takes the same form as (4.4), but with scale parameter = \/n/n2] 

(4.5) Y** ^ f{y**-y\T2). 

For the normal example, (4.5) is equivalent to generating 

Xl\...,Xll^Np{x\Ip) 

for given x* , and using the transformed variable Y** = ^/nX** . The two-step 
bootstrap probability with a pair of scales (ti,T2) is then defined by 

a2(2/,Tl,r2)=Pr{y**G7^;y,Tl,r2} 

ai{y* ,T2)f{y*;y,Ti)dy* , 
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where the integration is taken over the range of the components. We can 
write ai{y,Ti) = a2{y,Ti,0), because the conditional density of Y** con- 
verges to the point mass at y* by taking the hmit T2 — > 0. The two-step 
bootstrap might look similar to the double bootstrap of Hall (1992), but 
they are very different. We should generate thousands of Y** for given y* in 
the double bootstrap, but only one Y* in the two-step bootstrap. 
Let us consider two n2 values, 

(4.6) n2 = 6,15, 

for the normal example with parameter values (2.2). The two-step bootstrap 
probabilities are, for example, 

My, = 0.0359, My, yi) = 0.0205. 

Of course, they give ai{y, ^J^) and di(y, y^), respectively, in (3.2), be- 
cause 



a2 (y, Ti,T2) = ai{y,^J rf + r| ) 
for (1.1). For the exponential example with parameter values (4.1), however, 

"2(y, \/f,]/f) = 0-3063, a2{y, V^, Vb' ) = 0.1866 

are different, though very slightly, from di {y, \/^) = 0.2990 and di {y, \/ ^ ) 



0.1875, respectively, in (4.2). The difference of d2{y, ti,T2) from di{y, y Tj^ + 
for (1.3) is explained by 



(4.7) z2{y,n,T2)- zi{y,\/T^ + Tl 



(r2+r|)5/2 



We will use "=" to indicate that equality holds up to 0{n ^/'^) terms with 
error of order 0{n~^). Formula (4.7) and a revised regression model 

V — 2ft^^ 

(4.8) zi{y,Ti) = + {di-a)Ti 

for (1.3) are consequences of a more general argument with third-order ac- 
curacy shown in Section 7. 

The key idea in the two-step multiscale bootstrap is to estimate a by look- 



ing at the difference of a2(y,Ti,r2) from di{y, y + t"! )■ Once we compute 
Oii{y,Ti) and d2{y,Ti,T2) for several values of {ti,T2) by the one-step and 
two-step multiscale bootstrap, we can estimate v, di and d by fitting (4.7) 
and (4.8) to the observed bootstrap probabilities. A second-order accurate 
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p- value, denoted 02(2/), is then computed by using the estimated geometric 
quantities in the z-value 

(4.9) Z2{y) = v-di+a{l-v'^). 

This expression is shown to be equivalent to (2.3) up to 0(n~^/^) terms by 
using (4.8); zo{y) = v + di — a(l + 2-u^) and ZQ{fj{y)) = di — d. In the next 
section we will describe a procedure based on the above idea, as well as its 
refined version with third-order accuracy. 

It follows from (4.8) that the one-step multiscale bootstrap estimates v — 
2dv^ and di — d for the coefficients i) and c, respectively, under (1.3). Thus, 
zi{y) =v — di + 0(1 — 2v'^) = Z2{y) — dv"^, as well as %(y) = Z2{y) + 2di — 2d — 
dv'^, is first-order accurate in general. Since the difference ^2(2/) — zi{y) = dv'^ 
does not involve di, the one-step method adjusts the bias resulting from the 
curvature even if the normal model is misspecified. 

5. Three-step multiscale bootstrap resampling. We may repeat "step- 
ping" to obtain multistep-multiscale bootstrap probabilities so that we might 
be able to compute higher-order accurate p-values. This is the case, in fact, 
for going one step further, although the results are not known for yet further 
stepping. We introduce the three-step multiscale bootstrap for computing a 
third-order accurate p-value, denoted as{y), under (1.3) and (1.4). In the 
following argument, we first describe the procedure to compute a2iy), which 
helps understand that for as{y). 

The expression for Z2{y,Ti,T2) is obtained from (4.7) by substituting 

\Jti +T2 for Ti in (4.8). This is also expressed as 

(5.1) 52(y,Ti,r2) = (2(71, 72,73, Tl,r2), 
where the function ^2 on the right-hand side is defined by 

(5.2) C2(7i,72,73,n,r2) = 5171(1 + ^273) - 21lt£^. 

si7i 

Here si = (rf -|- tI )~^^^ and S2 = tI sf are functions of the scales, and the 
7i's are specified as functions of y under (1.3) and (1.4); 

(5.3) ^i = v — 2diP', 72 = t)(a — di), 73 = t)a. 
These 7i's are also used to express 

(5.4) z2(y)=7i(l + 73) + ?^, 

71 

which is equivalent to (4.9) up to 0{n~^/'^) terms. We calculate 02(2/, ''"i, ^2) 
for several values of (ti, r2) by the two-step multiscale bootstrap resampling. 
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and fitting tlie observed Z2{y,Ti,T2) = —^~^{a2{y,Ti,T2)) to the nonlinear 
regression model (5.1). Then the estimated 7j's are used to compute 02(2/) = 
^{-Z2{y)) from (5.4). 

This procedure is generalized for the three-step multiscale bootstrap re- 
sampling. A third-step replicate Y*** is generated for each y** by 

~/(y ,y ,T3) 

using the scale t^, and the three-step bootstrap probability is defined by 

a3{y,Ti,T2,T3) = Pr{y*** £ TZ]y,Tl,T2,T3} 

= j a2{y*,T2,T3)f{y*;y,Ti)dy*. 

Then, observed %(y, n, r2, ra) = —^~^{a3{y,Ti,T2,T3)) for several values of 
(''"i,T"2,T3) are fitted to the nonlinear regression model (3, defined by 

C3(7l> 72, 73, 74, 75, 76, Tl, T2, T3) 

(5.5) =71^1(1 +73'52 + 47IS2 -h75S3 +76S4) 

- (7i'5i)~H72 + 73'52 + 77|s2 + 74S2 + 375S3 + 376S4), 

where si, . . . , S4 are given by 

.1 = (rf + r| + r|)-V2, = (,2,2 + ^2^2 + ^2^2)^4^ 

^3 = (rMr| + riri + r,^(r| + r|)).?, .4 = (rMr|).?. 

The least squares estimates for the six 7j's are denoted by 71,..., 76. We 
then compute 0:3(2/) = <!>(— 23(1/)) by using the estimated 7i's in 

(5.6) Z3{y) = 71(1 + 73 + 47I + 76) + 7f 1 (72 + 7I/2 + 74 + 75)- 

Section 7 is mostly devoted to proving the third-order accuracy of 03(2/). 
The justification for the second-order accuracy of 02(2/) then immediately 
follows by ignoring 0{n~^) terms. As seen in (5.3), 71 is 0(1), and 72 and 73 
are 0(n^^/^). The rest of the three 0{n~^) geometric quantities are defined 
in Section 7.8. We do not have to know, however, the expressions of 7j's 
for computing 63(2/), because their values are estimated from the nonlinear 
regression, and the estimation error is only 0(n~^/^). 

It should be noted that there are other asymptotically equivalent expres- 
sions for ^3 and Z3 as functions of coefficients transformed from the six 7j's; 
we have shown the two different expressions for (^2 and Z2 as functions of 
either 71,72,73 or v,di,a. The expressions (5.5) and (5.6) are obtained by 
seeking simple ones. 
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0.145 + 0.008 - 0.018 - 0.0004 1 

+ } = 0.0509, 

1.328 



6. Examples. The two procedures in the previous section are apphed 
to the exponential example with parameter values (4.1). By the two-step 
multiscale bootstrap, the least squares estimates of 7j's are 

71 = 1.328 , 72 = 0.144, 73 = 0.137, 

and the corrected p- value is computed as 

a2{y) = 1 - <5{1.328(1 + 0.137) + ^} = 0.0528, 

which comes closer to the exact p- value aooiv) = 0.05 than ai{y) = 0.0753 
computed in (4.3). By the three-step multiscale bootstrap, the least squares 
estimates of the 7i's are 

71 = 1.328, 72 = 0.145 , 73 = 0.127, 

^4 = -0.018 , 75 = -0.0004, 76 = -0.036, 

and the corrected p- value is 

a3(y) = 1 - $|l.328(l + 0.127 + 0.065 - 0.036) 

4^ 

which is even better than a2{y) = 0.0528. 

In Table 2 p- values are computed for several parameter settings. The 
bootstrap probabilities are computed numerically {B = oo), but the stan- 
dard errors due to the bootstrap resampling are shown for B = 10,000. The 
first row corresponds to the normal model with (2.2), and the fourth row 
corresponds to the exponential model with (4.1). The following two rows for 
each are obtained by changing n = 10 to 100 and 1000. Similarly, the last 
six rows are obtained by changing ckcxd = 0.05 to 0.95. We observe that all 
the p-values tend to converge to cioo as n grows, and the corrected p-values 
are faster for convergence than oiQ. 

C(3iy,Ti,T2,T3) is Computed for all the combinations of (ri,r2,T3) values, 
as noted in the table; five (ti, 0, 0)'s, ten (ti, r2, 0)'s, and twenty (ti, r2, r3)'s. 
Therefore, the numbers of bootstrap probabilities are 5, 15 and 35, respec- 
tively, for ai{y), 02(2/) and a3{y). The nonlinear regression models are fitted 
to these bootstrap probabilities, and the least squares estimates of the ge- 
ometric quantities are calculated; each residual term is weighted inversely 
proportional to the estimated variance. For stable estimation, ridge regres- 
sion is also used; a penalty term J2i=iOJi'jf with small uji values is added to 
the residual sum of squares for minimization. 

For the exponential distribution, is kth. order accurate {k = 1,2,3), 
and, in fact, \ak — 0(qo\ becomes smaller as k increases in the table. It turns 
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Table 2 

p-values in percent (standard error) for the examples* 



Ridge regression 



n 


do 


Aabc 






d2 


ds 








Normal distribution {oaa = 5.00) 








10 


0.85 


7.75 


5.29 (0.61) 5.85 (1.81) 


7.03 (8.04) 


5.67 


(1.03) 


6.04 (1.13) 


100 


2.73 


5.25 


5.01 (0.37) 5.05 (1.16) 


5.08 (2.93) 


5.04 


(0.78) 


5.06 (0.97) 


1000 


4.12 


5.03 


5.00 (0.32) 5.00 (1.05) 


5.00 (2.22) 


5.00 


(0.72) 


5.00 (0.89) 








Exponential distribution (Aoo = 5.00) 








10 


11.15 


5.00 


7.53 (0.31) 5.28 (0.77) 


5.09 (0.95) 


5.77 


(0.60) 


5.13 (0.68) 


100 


6.73 


5.00 


5.90 (0.30) 5.03 (0.94) 


5.01 (1.50) 


5.25 


(0.67) 


5.04 (0.81) 


1000 


5.52 


5.00 


5.29 (0.30) 5.00 (0.98) 


5.00 (1.82) 


5.08 


(0.69) 


5.01 (0.80) 








Normal distribution ( 


doo = 95.00) 








10 


67.84 


92.33 


95.26 (0.18) 95.20 (0.41) 


95.02 (0.51) 


95.21 


(0.34) 


95.07 (0.37) 


100 


90.65 


94.74 


95.02 (0.24) 95.07 (0.84) 


95.09 (1.28) 


95.06 


(0.60) 


95.07 (0.70) 


1000 


93.91 


94.97 


95.00 (0.28) 95.00 (0.95) 


95.00 (1.72) 


95.00 


(0.67) 


95.00 (0.81) 








Exponential distribution (doo = 95.00) 






10 


98.78 


95.00 


97.99 (0.24) 94.48 (1.31) 


96.12 (7.39) 


95.60 


(0.81) 


96.48 (0.56) 


100 


96.49 


95.00 


95.95 (0.28) 94.97 (1.06) 


95.01 (2.71) 


95.24 


(0.72) 


95.14 (0.82) 


1000 


95.50 


95.00 


95.30 (0.29) 95.00 (1.02) 


95.00 (2.19) 


95.08 


(0.70) 


95.02 (0.81) 



*Tlic bootstrap calculation is replaced by integration numerically, and, hence, the number of 
bootstrap replicates is regarded as _B = oo. The standard errors in parentheses are calculated 
for the case of _B = 10'' by the local linearization of the nonlinear regression [Draper and 
Smith (1998)]. AU the combinations of r? e {f , f , 12, if, 12}, ^| g {f , if}, r| G {f , if} 
are used for the scales. The total numbers of bootstrap replicates are 55, 15-B and 35B, 
respectively, for di, d2 and da. For the ridge regression, the penalty weights are uji — iL>2 = 
and ujs = ■ ■ ■ — iJe ~ 0.01. 

out that |aabc — Oioo\ is almost zero here, because ciabc happens to be third- 
order accurate for the one-dimensional exponential distribution, as shown 
in Section 7.7. 

For the normal distribution, ai, 0,2 and 03 are third-order accurate, be- 
cause 73 = • • • = 76 = under (1.1), as shown in Section 7.8. This may explain 
why Idfc — ciool becomes larger as k increases in some of the rows. These four 
geometric quantities of zero value are estimated from slight differences of 
bootstrap probabilities, leading to unstable estimation as seen in the large 
standard errors. This is alleviated by ridge regression; even the worst case 
in the table = 6.04 it 1.13 may be allowed in practice. However, the total 
number of replicates is 350,000 for 03, almost comparable to that of the 
double bootstrap for achieving the same degree of the standard error. 

Although Qi is first-order accurate for (1.3), it is reasonably accurate 
even for the exponential model in the table. The total number of replicates 
is 50,000, yet the standard error is considerably smaller than that of 03. 
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Similar observation holds for the second-order accurate a^- The one-step, as 
well as two-step, multiscale bootstrap may provide a compromise between 
the number of replicates and the accuracy in practice. 

7. Asymptotic analysis of the bootstrap methods. 

7.1. A unified approach. Our approach to assessing the bootstrap meth- 
ods is not very elegant but rather elementary and brute-force. We explicitly 
specify a curved coordinate system along dTZ, which is convenient to work 
on the bootstrap methods. The density function of Y with respect to the 
curved coordinates is first defined for r = 1 in Section 7.2 and extended 
for r > in Section 7.3. We define a modified signed distance by altering v 
slightly, and its distribution function is given in Section 7.4. 

It turns out that the z-values of the bootstrap probabilities are special 
cases of the modified signed distance, and our approach gives an asymptotic 
analysis of the bootstrap methods in a systematic way. Using the result of 
Section 7.4, a third-order accurate pivot statistic is defined in Section 7.5, 
and the distribution functions of the bootstrap z-values are shown in Sec- 
tions 7.6 to 7.8, proving the main results of Section 5. 

The proofs of lemmas are given in Shimodaira (2004). We have used the 
computer software Mathematica for straightforward and tedious symbolic 
calculations; the program file is available from the author upon request. 

7.2. Tube- coordinates. In our curved coordinate system, a point rj is 
specified by two parts, a point on dTZ and the signed distance from it. This 
is an instance of the coordinate system used for the Weyl tube formula, 
and we call it tube-coordinates. Below we will define the coordinate system 
explicitly, and show the expression of the density function of Y in terms of 
the tube-coordinates. We take an approach similar to that of Kuriki and 
Takemura (2000). 

The density function of the exponential family of distributions is expressed 

as 



where 9 = {9^, . . . , 9^) is the natural parameter vector. We denote (7.1) by 
f{y;r]) using the expectation parameter vector rj = (rji, . . . ,rjp) = E{Y), the 
expected value of Y. The change of variables <-> is one-to-one, and is given 
by rji = dip/d9'^, 9^ = dcp/drji, i = 1, . . . ,p, where the potential function (j){ri) 
is defined from the cumulant function ^(0) by (j){r]) = maxg{9^r]i — ^(0)}. 
The metric at t] is denoted as 





d-rji dr]j 
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and the derivatives of (p at rj = are denoted as 



d(p{rj) 



drji 



d^(j){ri) 



drji drjj 



d^(t){ri) 



'JiH ^'Ij 



drji drjj drjk 



and so on. 



Since the exponential family is not uniquely expressed up to affine transfor- 
mation, we assume without loss of generality that = and (p''^ = 5ij, where 
6ij takes value one when i = j, otherwise zero. In other words, E(Y) = and 
cov(y), the covariance matrix of Y , is Ip at 9 = 0. We make our asymptotic 
argument local in a neighborhood of r/ = by assuming the local alternatives. 
The smooth surface dTZ of the region TZ is specified locally around r/ = 

by 

r]a{u)=Ua, a = l,...,p-l; r]p{u) ^ -d^^UaUb - e'^^'^UaUbUc, 

where u = (ui, . . . , Up-i) is the (p — l)-dimensional parameter vector to spec- 
ify a point T]{u) on dTZ. TZ is specified locally by r]p < rjp{u). It follows 
from the argument below equation (2.12) of Efron and Tibshirani (1998) 
that dp'' = 0(n-i/2) ^nd e'''"= = 0(n-^), and similarly, (p^^^ = 0(n-^/2) 

Let Bf{u) = drji/dua-, i = I, . . . ,p, be the components of a tangent vector 
of the surface for a = 1, . . . ,p — 1. They are given explicitly as 

B^{u) = 6ab, 6 = l,...,p-l; S^(n) « -2(i"^b - Se'^^^feUe, 

and the metric in the tangent space is given by 

(7.2) « 6ab + 'f'W 

+ {W'^d^'^ - 2d''^(t)^'^^ - 2d^'^(t)''^P - d'"^(P''''P + ^(l)^^'"^}ucUd, 

where (j)'^{r]{u)) 6ij + (p'^'^Ua + {-d'^'^cp'^P + ^(p'^'^^^juaUb. Let Bf{u), i = 
1, . . . ,p, be the components of the unit length normal vector orthogonal to 
the tangent vectors with respect to the metric such that 

cP^^{v{u))Bt{u)BP{u) =0, a = 1, . . . ,p - 1; 

cl>^^{rj{u))Bf{u)BP{u) = l. 

The components are calculated explicitly as BP{u) ~ {2d"''' — (p'^''P)ub + {'ie°''"^ + 

^ab^cpp dbc^app _ 2d^d^acd _^ ^abd^cdp _^ l^abp^cpp _ l^abcp^^^^^^ 



U 



1 _ IfP'^PPua + {-2d"'=d^'= + ld"''(j>PPP + + ^(P"PP(P''PP - \(P"''PP}UaUb. 

Let u be a scalar, and {u,v) be a p-dimensional vector. We consider repa- 
rameterization defined by 

(7.3) r]i{u,v) =r]i{u) + Bf{u)v, i = l,...,p, 
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and assume rj {u,v) is one-to-one at least locally around r] = 0. {u,v) 
gives the tube-coordinates of the point rj. The boundary dTZ is expressed 
simply by ?; = 0, and the region 7^ is u < 0. {u,v) is used for indicating the 
parameter value tj = r]{u,v), or the observation y = ri{u,v). When there is a 
possibility of confusion, we may write y ^ {u,v) instead of r] ^ {u,v). 

Since the normal vector is orthogonal to the surface, rj{u) = rj{u, 0) G dTZ 
is the projection of rj{u,v) onto dTZ; u is the maximum likelihood estimate 
under the restricted model specified by dTZ. r]{u, 0) is denoted by f}(?/) in 
Section 1 as a function of y. v is the signed distance mentioned for (1.1) in 
Section 3. 

V is also related to the signed likelihood ratio R [McCullagli (1984) and 
Severini (2000)] hy R^v + ^^pPv^ + {^^ppp - ^(^PP)2}^)^ where ^pp 
and (ffPPP are the third and fourth derivatives to the normal direction eval- 
uated at rj{u,0), instead of ry = 0. This third derivative is associated with 
the acceleration constant. For the acceleration constant a, the formula a = 
— ^(f>PPP is obtained directly from equation (2.9) of DiCiccio and Efron (1992), 
or by using equation (6.7) of Efron (1987) and d^ip/dO' 89^ 89^ = -cp'^''. 
The expression for the density function of iU,V) is obtained from f{y;r]) 
by change of variables, as shown in the following lemma. 

Lemma 1. Let Y ~ f{y; rf) he the exponential family of distributions with 
r] = E{Y). Without loss of generality we may assume that cov(y) = Ip at 
rj = and that the true parameter value is specified by r] = {0, . . . ,0, X) for 
some X, that is, r]a = 0, a = 1, . . . ,p — 1, rjp = X, or, equivalently, u = 0, 
V = X using the tube- coordinates {u, v) r]. Let f{u, v; X) be the joint density 
function of {U,V) Then, ignoring the error of 0{n~^/'^), we obtain 

log /(u, v; X) w g{v, X) + g"-{v, X)ua + g^^iv, X)uaUb 

(7.4) 

+ g°''"'{i), X)UaUbUc + g''^'"^{v, X)UaUbUcUd, 

where the five functions on the right-hand side are defined by g{v,X) = 
-iplog(27r) - i(t) - A)2 - + i((/'^^'^)2 - ^(j^PPX^ - ^(jfPPPX^ + {2d'"' - 

l^aap l^pp Ifj^PPPX'^ + ^(f)PPPPX3jy + {_2(d«'')2 + 2d''^ (p^'^P - f ((^''^P)^ - 
i((/."PP)2 - l{(j)PPPf + \(PPPPP + 1(/,"«'PP}{;2 _ ^(j)PPPv3 - ^(jyPPPPi,^^ g^[y^ \) = 
1 (pabb 1 ^appy2 1 ^apppx^ { - A - d"'' (j)'"'"' + hd''^(j)''PP + (p'^PPd'"'' - Icj)'"'"' d''" + 

l^abp^bcc _ S^abp^bpp _^ l^app^bbp _ 3^app^pp _^ l^abc^bcp _ l^abbp ^ l^appp _^ 
Q^abb _^ fiab^bppx2 _ l^abp^bppx2 _ 1 ^app^pp;y2|^ _^ {_(iab^bpp ^ ^^abp^bpp _^ 

l0appp}^3^ 5«f'({;,A) = -\5ab- d''^X-\d''^(t)''''P + \(t)"-^''^ -\(t)''"^4>^"^ + 
2(pc^bc _ 2cfac^fecp _ l^ab^ppx2 + + 1 ^"fep _ (^2d'"'d^'' - id'^*</.PPP+ - 

^^acp^bcp _ 3^<ipp0f'PP)A}{), 5'^'"=({), A) = -^ct)"'"' - e^^^X + {-le"^" + ^cp^^^P - 

3^ab^cppj^^ad^bcd_l^abd^cdp_l^abp^cpp^^^ gabcd(^^^X) = -^d""^ d""^ + \(j)''^P d'"^ - 
1 jjubcd 

24:9 
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7.3. Changing the scale. We define a density function f{y',r],T) witli 
mean 77 and scale r > by modifying /(y; rj). Here r is regarded as a known 
constant, whereas r/ is a unknown parameter vector. Let (j){rj,T) be the po- 
tential function of f{y]rj,T), and (f){r]) be that for f{y;ri). Since the density 
function is defined by specifying the potential function, the following equa- 
tion gives a definition of /(y;r/,r): 

(7.5) 4>{v,T) = ^{r])/T'. 

This /(y;r/,r) comes naturally from the multiscale bootstrap resampling. In 
fact, the potential function of the replicate Y* is (j){r],T) = ||77|p/(2r^) for 
the normal example (2.1) of Section 2, and that is (j){r], r) = —n{l + logr]) /r^ 
for the exponential example of Section 4, and thus both agree with (7.5). 
The same applies to the exponential family, in general, as shown below. 

Lemma 2. Let X be a p- dimensional random vector of the exponential 
family. We assume that Y is expressed as a sum of m independent s such 
that Y = ^/n{Xi + •••-!- Xm)/m for m > 0, and that the density function is 
f{y; rj) when m = n. Then Y ~ /(y; r], r) with t = yjnjm for r > 0. 

We continue to use the tube-coordinates defined by the reparameteri- 
zation 77 •s-^ [u^v) of (7.3). By altering the potential ^{r]^V) to (j){r],T), the 
metric, as well as the tube-coordinates, should have changed if we go back 
to the specification of 'r](u) and BP(u) given in the previous section. How- 
ever, we continue to use the specification with r = 1 for any r > 0, so that 
the reparameterization rj <-> {u,v) does not depend on r. 

Lemma 3. Let f{u,v;X) be the joint density function of {U,V) ^ Y 
given in Lemma 1, and f{u,v;X,T) be that corresponding to f{y','r],T) with 
scale r > 0. Then the expression of log f{u,v;X,T) is obtained from (7.4) by 
changing {u,v) to 

(7.6) u = u/t, v = v/t, 

by adding the logarithm of the Jacobian log(l/T^) to (7.4), and replacing 
(ffjk ^ d'^^ , e'^^'^ and X, respectively, with 

(7.7) 

^ ^^ab^ ~abc ^ ^2^abc^ ~^ ^ ^^^^ 

7.4. Modified signed distance. We consider yet another transformation of 
the coordinates for expressing the bootstrap z-values in modified v values. 
Let w he a scalar variable defined formally by the series 

00 00 
(7.8) U) = V + '^CrV'' + Uc'^b'}.v'', 

r=0 r=0 
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where denotes the rth power. The coefficients are Cr 
bj. = 0{n~^), and their expressions are specified later. We assume the trans- 
formation {u,v) <-> {u,w) is one-to-one at least locally around {u,v) = 0. By 
inverting the series in (7.8), we also have 

oo oo 

(7.9) V = w — Crw'^ — Uc b^W^ , 

r=0 r=Q 

where Cr = Cr — J2l=oi^~ ■^ + ^)'^r-s+iCs, and b^ = b^- The coefficients are Cr = 
0(n~^/^) and b^ = 0{n~^). Let W be the random variable corresponding to 
w; the observed value w is defined by (7.8) but using the observed {u,v) 
instead of (n, v). 

We call w a modified signed distance characterized by the coefficients 
b'^, Cr] w reduces to v when all these coefficients are zero. The z- values of 
the bootstrap probabilities are represented as w by appropriately specify- 
ing the coefficients. The following lemma plays a key role in studying the 
distributional properties of the bootstrap probabilities. 

Lemma 4. Let us assume that the distribution ofY in the tube- coordinates 
is specified by {U,V) ^ f{u,v; \,t), and the coefficients in (7.9) are of or- 
der 6^ = O(n-i) for r > 0, cq = 0{n-^/^), d = 0{n"^), C2 = 0{n~^/^), 
C3 = 0{n~^) and Cr = 0{n~^/^) for r > 4. We define Zc{w;X,t) from the 
distribution function of the modified signed distance W as 

Pv{W <w} = ^{zc{w; A, r)). 

Then the Zc-formula is, ignoring the error of 0{n~^^^) , expressed as 

(7.10) Zc{w;X,T)^T~'^g^{w,X) -^Tg+{w,X), 

where g- (w, X) = {w- X) -cq- ic/^^PfA^ + I^ppXw + {l(tfPP- C2)w'^ - IcQcffPPX - 
{ci + lco(i)PPP}w + {l{(l)^PP)^ + ^{(l)PPP)^-l(t>PPPP}X^ + {-l{r^^ 

{_^(0PPP)2 + l_^PPPP _ lc2(/>PPP}Au;2 + {-^{^PPP)^ + ^(f>PPPP - \C2(FPP - 

c^}w^, and g+{w, A) = -{d'"' + ^^pp) + {(d"*)^ - d'^'^^'^^P + ld'"'(l)PPP + ^{^''''p)^ + 

1 (<^app)2 ^ 13 (^PP)2 _ l^aapp _ l^PPP}^'; + {((i»'')2 _ ^^^^qyPPP + ^cP^PPf + 

^{(t>PPP)^ — ^(t)PPPP}X. Note that the Zc-formula does not involve the coef- 
ficients b'^., and that the distribution function of W is characterized by the 
coefficients Cr with third-order accuracy. The index c of Zc indicates the co- 
efficients Cr- 

The true parameter value is assumed to be (0, A) in the (u, i))-coordinates 
for (7.4) and (7.10). If we alter the true parameter value to arbitrary (u,t;) 
with ti 7^ 0, the expression changes as well, and <I>~^(Pr{Ty < li)}) is denoted 
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as Zc{w; u, v, r), which reduces to Zc{w; 0, A, r) = Zc^w; A, r) when u = and 
v = X. 

Zc{w;u,v,t) is used for representing the bootstrap probabihties in par- 
ticular. The simple bootstrap probability is, for example, ao{y) = Fr{V* < 
0;y} = ^{zc{0;u,v,l)) with all Cr = 0. The expression of Zc{'w*;u,v,t) is 
obtained from (7.10) by changing the origin to 'r]{u). 

Lemma 5. Let Y* be a replicate ofY distributed conditionally as Y* ~ 
f{y*',y,T) with mean y and scale t, and W* be the corresponding modified 
signed distance. Let us denote the conditional distribution of W* given y as 
Ft{W* < w*-,y} = ^{zc{w*-,u,v,t)). Then the expression of Zc{w*]il,v,T) 
is obtained from (7.10) by replacing w, \, (p^^P and di = d"", respectively, 
with w* , V, 

(7.11) (^PP w (jfPP + {3(p^PP{2d!"' - (t)^^P) - ^(P^PPcpPPP + (P^PPPjuc and 

(7.12) di « d"" + - 'i^V"^" + 3e'"'^}uc. 

Note that 0{n~^) terms change only 0(n~^/^). For example, d2 = {d""^)"^ 
would be replaced with d2 , but d2~d2. 

7.5. Pivot statistic. Although the exactly unbiased p- value may not exist 
in general, a third-order accurate p-value can be derived under (1.3) and 
(1.4). Let Y* ~ f{y*]fl{y)^ 1) be a replicate generated with mean f/(y) instead 
of y, and a.oo{y) be defined as the probability of the corresponding signed 
distance V* being greater than or equal to the observed value v] 

aoo(2/)=Pr{F*>?);7}(y)}. 

This is the exact p-value for the normal example of Section 2 and for the 
exponential example of Section 4. We will show that aoo(y) is, in fact, third- 
order accurate under (1.3) and (1.4). 

First, Zoo{y) = —^~^{oioo{y)) is expressed by the Zc-formula of Lemma 5. 
From the definition, ioo(y) = Zc{v] u, 0, 1) with all = and, thus, 

Zoo{y)^v-{di + l4>'n + h^''''v' 
(7.13) 

+ ^{^"''p)'^ + i((/.'^pp)2 + m(j)ppp)^ - i^'^'pp - ici)pppp}v 

By comparing (7.13) with (7.8), we find that Zoo{y) can be expressed as 
w with coefficients cq = -d"" - l(j>PPP, ci = {d^'^f - d''^(l)''^P + \d'"'(tPPP + 
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l(jfppv^ Then the distribution function of ioo(y) is obtained immediately 
from Lemma 4 as shown below. 

Lemma 6. Let us consider a statistic 

Zq{y) ~ Zoo{y) + qo + qiv + q2v'^ + qsv^ + Ucg'^iv), 

where the coefficients are go = qi = 0{n~^), q2 = 0{n~^^'^) and 

(73 = 0{n~^), and g'^{v) = 0{7i~^), c = 1,. . . ,p — 1, representing arbitrary 
polynomials of v. The index q of Zq indicates the coefficients. Assuming 
([/, y) ~ /(n, -C; A, 1), the distribution function of Zq{y) is expressed as 

Vl{Zq{Y)<X-A} 

« «>[X - A - go - \4^^^\^ + l(t)P^^\x - 92X2 

(7.14) 

+ {-gi - 2q2(d^'' + i<APPP - go)}x + {-\{4>'''"'f + ^^PPPPjA^rE 
+ {\^'""'q2 + 2gi - q^}x' + {i(<A'^^^)2 + U^''^? " W""}^'' 

+ {-^('Z'^^^)' + - \(\?^^q2\\A- 

For A = 0, the distribution function is Pr{zg(y) < x;0} ^\x — qo — 
q2X^ + {-qi - 2q2{d'"' + ^c/^pp - qo)}x + {i(/^PPg2 + - gajx^]. In par- 
ticular, Pv{zoo{Y) < x;0} ^{x) and, thus, Zoo{y) is a third-order accurate 
pivot statistic. We obtain Pr{Q;oo(^^) < a;rj} ^ a for G 37^, proving the 
third-order accuracy of aoo(y)- 

The reverse of the above statement also holds. aq{y) = ^{—Zq{y)) is a 
third-order accurate p- value if and only if go ~ gi ~ §2 ~ <?3 ~ 0. If we confine 
our attention to a-q{y) defined only from v and the geometric quantities 
^ab^ ^abc^ ^ij ^ ^ijk ^zjki evaluated at ?7(y), then Ucg^{v) in Zq{y) comes 
only from g^'s by the replacements shown in Lemma 5. Thus, aq{y) is a 
third-order accurate p- value if and only if aq{y) ~ ctooiu)- Similarly, aq{y) 
is second-order accurate if and only if go = 92 = and, thus, aq{y) = aoo(y)- 

-Soo(y) is equivalent to other pivots in the literature up to 0{n~^) terms. 
Under (1.1) and (1.4), 0*^'= = (j)'^''^ = and, thus, (7.13) reduces to Zoo(y) » 
V — di + d2V, giving (3.8), the pivot of of Efron (1985). Under (1.3), the 
modified signed likelihood ratio [Barndorff-Nielsen (1986) and Barndorff- 
Nielsen and Cox (1994)] has been known as a third-order accurate pivot, 
and it is expressed as R* = R + {1 / R) log{U / R) in the notation of Severini 
[(2000), page 251], where U is defined using the log-likelihood derivatives. A 
straightforward calculation shows that U k,v — div"^ + {^{d"''^)'^ + d'^'^d"-^ — 

yaapp _ ^ab^abp ^ l(^ai,p)2 ^ l(^app)2 ^ 1(<^PP)2 _ _^^PPP}^)3^ ^ud that 

R* 

~ Zooiu) ill tliG moderate deviation region. 
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7.6. Accuracy of the bootstrap probability. Since the event Y* £ TZ is 
equivalent to the event V* < 0, the z-value of the bootstrap probabihty with 
scale r is expressed by the Zc-formula of Lemma 5; zi{y,T) = —Zc{0;u,v,t) 
with all Cr = 0. From (7.10), we obtain a refined version of (4.8), erring only 

0(n-3/2), 

zi{y,r) ^ T-'[v + \4PPPv^ - {i(<A'^^^)2 + U^n"" - W^nA 

(7.15) +r[{d, + \^^PP) 

- {{d''^? - Id^^'clfPP + |(<^''^^)^ + ^(<^PP)^ - ^^PPP}v]. 

It follows from (7.15) that rzi {y, r) is expressed as w and, thus, tzi {y, r) ^ 
Zq{y) by choosing the coefficients appropriately. They are cq = {d°'°' + |(^^^^)r2, 

and C3 = -1(0»PP)2 _ A(<^PP)2 + i^PPP for or, equivalently, qo = {I + 
T^)id'"' + l(j)PPP), qi = -(1 +r2)(d'^*)2 + d^^^^^P + IcP'^'PP - i((/>'^''P)2 - 1(4 + 
r2)(</)'^PP)2+ i(-l + r2)d'^»<^PP--L(l3 + 5r2)(0PPP)2 + -L(3 + r2)<^PPP, 53 = 
^(j)PPP^ 53 = _ 1 ((^'^PP)2 _ -L(0PPP)2 + ^^pppp for zg(y). The distribution func- 
tion of Tz{y,T) is obtained from (7.10) or (7.14). In particular, the distribu- 
tion function of zo{y) = zi{y, 1) under A = 0, r = 1 is 

^<^[x- {2d'''' + ^cfPP) - \(t)PPPx^ 

(7.16) -t- {2(d'^*)2 - (i'^V"''^ + \d''''(lfPP + i((/.°*P)2 

+ |(^"PP)2 + li(,/,PPP)2 _ 1</,««PP _ 

showing the first-order accuracy of ao(y). 

Remark A of Efron and Tibshirani (1998) discusses a calibrated boot- 
strap probability, denoted adoubic(y) here, using the double bootstrap of 
Hall (1992). Similarly to the two-level bootstrap, thousands of Y* are gen- 
erated around r/(y). Then do(2/*) is computed for each y* . The expression 
of 5doubic(y) = $"MPr{%(r*) < zo{y);fi{y)}] is obtained from (7.16) by the 
replacements of Lemma 5, and a straightforward calculation shows that 
^double (y) ~ 5oo(y), proving the third-order accuracy of Q;double(y)- 

7.7. Accuracy of the two-level bootstrap. The expression of ZQ{y) is ob- 
tained from (7.15) by letting r = 1, and zo(^(y)) ~ di + \(j)PPP is obtained 
from it by letting {) = 0. By substituting these expressions, as well as a = 
— ^(j)PPP for those in (2.3), we find that fabc(y) is expressed as w, or, equiva- 
lently, Zg{y) with coefficients go = 92 = 0, qi = -2{d^^f + ^(jy^'^PP + - 
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J-fj^ppp^ The distribution function is then obtained from Lemma 6. For A = 0, 
it becomes 

(7.17) Pr{iabc(^) <x;0}^ $(x - qix - q^x^), 

showing the second-order accuracy of aabc(y)- 

For the exponential example of Section A, p=l, <p^^^ = — 2/-^/n, (p^^^^ = 
6/n and all the other quantities in qi and q^ are zero. Therefore, qi = q3 = 
0, and Zshciu) turns out to be third-order accurate, explaining the high 
accuracy of aabc(y) observed in Table 2. 

7.8. Accuracy of the multistep-multiscale bootstrap. Using the expres- 
sions (7.4) and (7.15), the expression of Z2{y,Ti,T2) is obtained by the inte- 
gration 

(7.18) ~Z2{y,n,T2) = ^-^y ^~zi{y*,T2))f{y*;y,Ti)dy*y 

By repeating the same integration using Z2{y* ,T2,T3) instead of 2i(y*,r2), 
we obtain the expression of ^3(y, ti, r2, rg) as given below. 

Lemma 7. Let us define the following six geometric quantities using the 
derivatives evaluated atr] = 0:ji = X + ^X^^pp + \^{-li(l)^PP)^ - ^{(j)PPP)^ + 
l<j)PPPP}, 72 = A{-d'^" - l(J)PPP} + X^{{d<'^)^ - y'"'(l)PPP + l{(j)^PPf + j^{(j)PPP)'^ - 

±(j)PPPP}, 73 = -^X(l)PPP + X^{l{cP<^PPf + ^(cj)PPP)^-lcl)PPPP}, 74 = A2{-d«''<^«^P + 
^^aa^pp 1 ((^afep)2 _^ 1 (0app)2 2 (^^pp)2 _ ^^aapp _ 1 ^PPP} ^^^ = X^{-\ {(p'^PP)'^ - 

i(<^PP)2 + ^^PPPP} and 76 = A2{-|((/>"pp)2 - ^{(jfPP)^ + ^ctPPPP}. Those 
evaluated at 7]{y), denoted 71,..., 70, are obtained by replacing X, (jfPP and 
d"""" , respectively, with v, (7.11) and (7.12) as shown in Lemma 5. Then we 
have 

(7.19) ^3 (y , n , T2 , ra ) W Cs (71 , 72 , 73 , 74 , 75 , 76 , n , T2 , Ts) 

using the C,3-function of (5.5). Since (7.19) errs only 0{n~^/'^) for any val- 
ues 0/ (ti, r2, Tg), the nonlinear regression for three-step multiscale bootstrap 
probabilities in Section 5 estimates ^iS up to 0{n~'^) terms. 

If we define z-i{y) of (5.6) using the 7j's defined above, we can easily verify 

(7.20) h{y)-Zoo{y) 

by comparing (5.6) with (7.13). This proves the third-order accuracy of 63(2/) 
under (1.3) and (1.4). 

For the multivariate normal model of (1.1), </>(r/) = ||r7p/2 and, thus, 
^ijk _ ^ijki _ This implies 73 = • • • = 76 = 0, proving the third-order ac- 
curacy of a.i{y) and a2{y) under (1.1) and (1.4). 
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